function hfile = vmp_PerformFDR(hfile, opts)
% VMP::PerformFDR  - perform FDR correction (set properties)
%
% FORMAT:       [vmp = ] vmp.PerformFDR([opts])
%
% Input fields:
%
%       opts        1x1 struct with optional fields
%        .bonfval   alternative BonferroniValue (override map value)
%        .mapsel    map selection, default: all supported maps in VMP
%
% Output fields:
%
%       vmp         VMP object

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% settings
global bvqxfile_config;
fdrc = bvqxfile_config.settings.Statistics.FDR.Thresholds;

% check arguments
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmp')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if bc.NativeResolutionFile == 0
    return;
end
if nargin < 2 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'bonfval') || ...
   ~isa(opts.bonfval, 'double') || ...
    numel(opts.bonfval) ~= 1 || ...
    isinf(opts.bonfval) || ...
    isnan(opts.bonfval) || ...
    opts.bonfval < 2
    opts.bonfval = [];
else
    opts.bonfval = round(opts.bonfval);
end
if ~isfield(opts, 'mapsel') || ...
   ~isa(opts.mapsel, 'double') || ...
    any(isinf(opts.mapsel(:)) | isnan(opts.mapsel(:)) | ...
        opts.mapsel(:) ~= fix(opts.mapsel(:)) | opts.mapsel(:) < 1)
    opts.mapsel = 1:numel(bc.Map);
else
    opts.mapsel = unique(min(opts.mapsel(:), numel(bc.Map)));
end

% iterate over maps
for mc = 1:numel(opts.mapsel)
    
    % get map handle
    map = bc.Map(opts.mapsel(mc));
    
    % continue if map type is not supported
    if ~any([1, 2, 4] == map.Type)
        continue;
    end
    
    % get map values we care about
    if ~isempty(opts.bonfval)
        bonfval = opts.bonfval;
    else
        bonfval = map.BonferroniValue;
    end
    msk = true(size(map.VMPData));
    if mod(map.ShowPositiveNegativeFlag, 2) == 0
        msk(map.VMPData > 0) = false;
    end
    if map.ShowPositiveNegativeFlag < 2
        msk(map.VMPData < 0) = false;
    end
    mpv = map.VMPData(msk);
    if isempty(mpv)
        bc.Map(opts.mapsel(mc)).NrOfFDRThresholds = 0;
        bc.Map(opts.mapsel(mc)).FDRThresholds = zeros(0, 3);
        continue;
    end
    
    % computation depends on map type
    switch (map.Type)
        case {1} % t-map
            fdrv = [fdrc(:), applyfdr(abs(double(mpv)), 't', ...
                fdrc(:), map.DF1, [], true)];
        case {2} % r-map
            fdrv = [fdrc(:), applyfdr(abs(double(mpv)), 'r', ...
                fdrc(:), map.DF1, [], true)];
        case {4} % F-map
            fdrv = [fdrc(:), applyfdr(abs(double(mpv)), 'F', ...
                fdrc(:), map.DF1, map.DF2, true)];
    end
    
    % store
    bc.Map(opts.mapsel(mc)).NrOfFDRThresholds = size(fdrv, 1);
    bc.Map(opts.mapsel(mc)).FDRThresholds = fdrv;
end

% put content back in object
bvqxfile_setcont(hfile.L, bc);
